# Replication code for "Can Deliberative Minipublics Influence Public Opinion? Theory and Experimental Evidence"
# Goal: replicate the tables and figures in the main text
# Dependencies: dlexpdata.Rdata

library("survey")
library("xtable")
library("boot")
library("gplots")
library("Zelig")
library("ZeligChoice")

load("dlexpdata.Rdata")

#-------------------------- Table 1 --------------------------#

## Frequency table

# cross-tab
T1.freq.a <- table(dlexpdata$video, dlexpdata$party.cues)

# add row totals
T1.freq.b <- cbind(T1.freq.a, rowSums(T1.freq.a))               
               
# add column totals
T1.freq <- rbind(T1.freq.b, colSums(T1.freq.b))

## Percentage table

# cross-tab
T1.per.a <- prop.table(T1.freq.a) * 100

# add row totals
T1.per.b <- cbind(T1.per.a, rowSums(T1.per.a))               

# add column totals
T1.per <- rbind(T1.per.b, colSums(T1.per.b))

# print Table 1

colnames(T1.freq) <- colnames(T1.per) <- c("No partisan cue", "Partisan cue", "Total")
rownames(T1.freq) <- rownames(T1.per) <- c("No deliberative cue", "Deliberative cue", "Total")

xtable(T1.freq, digits = 0)
xtable(T1.per, digits = 1)

#-------------------------- Figure 2 --------------------------#

sdesign<-svydesign(id=~1, weights=~psweights.trim, data=dlexpdata)

dlexpdata$agree.num <- NA
dlexpdata$agree.num[dlexpdata$agree == "Approve strongly"] <- 5
dlexpdata$agree.num[dlexpdata$agree == "Approve somewhat"] <- 4
dlexpdata$agree.num[dlexpdata$agree == "Neither approve nor disapprove"] <- 3
dlexpdata$agree.num[dlexpdata$agree == "Disapprove somewhat"] <- 2
dlexpdata$agree.num[dlexpdata$agree == "Disapprove strongly"] <- 1

### First make tables used as input for Figure 2 and Table 2

## Condition: no partisan cues

dlexpdata$treat1vs2 <- NA
dlexpdata$treat1vs2[dlexpdata$video == 0 & dlexpdata$party.cues == 0] <- 0
dlexpdata$treat1vs2[dlexpdata$video == 1 & dlexpdata$party.cues == 0] <- 1

tplot1 <- prop.table(svytable(~dlexpdata$treat1vs2 + dlexpdata$agree.num, sdesign, exclude=NULL, na.action=na.pass), 1)[1:2,]
tplot1 <- rbind(tplot1, tplot1[2,] - tplot1[1,])
colnames(tplot1) <- rev(c("I'm not sure", "Approve strongly", "Approve somewhat", "Neither approve nor disapprove", "Disapprove somewhat", "Disapprove strongly"))

# obtain confidence intervals using bootstrapping

func.props.boot.1vs2 <- function(data, indices) {
  d <- data[indices,] # allows boot to select sample
  tplot1 <- prop.table(table(d$treat1vs2, d$agree.num, useNA = "always"), 1)[1:2,]
  tplot1 <- rbind(tplot1, tplot1[2,] - tplot1[1,])
  return(tplot1)
}

set.seed(1000)

props.boot.1vs2 <- boot(data = dlexpdata, statistic= func.props.boot.1vs2,R=1000, sim = "ordinary", simple = TRUE, weights = dlexpdata$psweights.trim)

tplot1.ci.low <- matrix(apply(props.boot.1vs2$t, 2, quantile, p = 0.025), nrow = 3, byrow = F)
tplot1.ci.high <- matrix(apply(props.boot.1vs2$t, 2, quantile, p = 0.975), nrow = 3, byrow = F)

# create table with standard errors

tplot1.se <- matrix(apply(props.boot.1vs2$t, 2, sd), nrow = 3, byrow = F)

tplot1.w.se <- t(rbind(tplot1, tplot1.se[3,])) * 100

colnames(tplot1.w.se) <- c("Condition I", "Condition II", "Difference (II-I)", "s.e.")
rownames(tplot1.w.se) <- c("Disapprove strongly", "Disapprove somewhat", "Neither approve nor disapprove", "Approve somewhat", "Approve strongly", "I'm not sure")

## Condition: partisan cues

dlexpdata$treat3vs4 <- NA
dlexpdata$treat3vs4[dlexpdata$video == 0 & dlexpdata$party.cues == 1] <- 0
dlexpdata$treat3vs4[dlexpdata$video == 1 & dlexpdata$party.cues == 1] <- 1

tplot2 <- prop.table(svytable(~dlexpdata$treat3vs4 + dlexpdata$agree.num, sdesign, exclude=NULL, na.action=na.pass), 1)[1:2,]
tplot2 <- rbind(tplot2, tplot2[2,] - tplot2[1,])
colnames(tplot2) <- rev(c("I'm not sure", "Approve strongly", "Approve somewhat", "Neither approve nor disapprove", "Disapprove somewhat", "Disapprove strongly"))

# obtain confidence intervals using bootstrapping

func.props.boot.3vs4 <- function(data, indices) {
  d <- data[indices,] # allows boot to select sample
  tplot2 <- prop.table(table(d$treat3vs4, d$agree.num, useNA = "always"), 1)[1:2,]
  tplot2 <- rbind(tplot2, tplot2[2,] - tplot2[1,])
  colnames(tplot2) <- rev(c("I'm not sure", "Approve strongly", "Approve somewhat", "Neither approve nor disapprove", "Disapprove somewhat", "Disapprove strongly"))
  return(tplot2)
}

set.seed(2000)

props.boot.3vs4 <- boot(data = dlexpdata, statistic= func.props.boot.3vs4,R=1000, sim = "ordinary", simple = TRUE, weights = dlexpdata$psweights.trim)

tplot2.ci.low <- matrix(apply(props.boot.3vs4$t, 2, quantile, p = 0.025), nrow = 3, byrow = F)
tplot2.ci.high <- matrix(apply(props.boot.3vs4$t, 2, quantile, p = 0.975), nrow = 3, byrow = F)

# create table with standard errors for Table 2

tplot2.se <- matrix(apply(props.boot.3vs4$t, 2, sd), nrow = 3, byrow = F)

tplot2.w.se <- t(rbind(tplot2, tplot2.se[3,])) * 100

colnames(tplot2.w.se) <- c("Condition III", "Condition IV", "Difference (IV-III)", "s.e.")
rownames(tplot2.w.se) <- c("Disapprove strongly", "Disapprove somewhat", "Neither approve nor disapprove", "Approve somewhat", "Approve strongly", "I'm not sure")

## create Figure 2

pos1 <- tplot1
pos1[,1:5] <- TRUE
pos1[,6] <- FALSE

cols1 <- matrix(c("grey26", "grey56")[pos1 + 1], nrow = 3, byrow = F)

pos2 <- tplot2
pos2[,1:5] <- TRUE
pos2[,6] <- FALSE

cols2 <- matrix(c("grey26", "grey56")[pos2 + 1], nrow = 3, byrow = F)

altnames <- colnames(tplot1)

par(mfrow=c(2,2)) 
par(oma=c(3,18,0,1),mar=c(1,2,6,1)) 

barplot2(tplot1[1,1:6], plot.ci = TRUE, ci.l = tplot1.ci.low[1,1:6], ci.u = tplot1.ci.high[1,1:6], xlim = c(0,0.5), horiz = TRUE, cex.names = 1.2, cex = 1.2, col = cols1[1,], las = 2, axisnames = TRUE)
mtext("Condition I", side = 3, line = 0, cex = 1.2)
mtext("No partisan cue", side = 3, line = 2, cex = 1.4)
mtext("No deliberative cue", side = 2, line = 16, cex = 1.4)

barplot2(tplot2[1,1:6], plot.ci = TRUE, ci.l = tplot2.ci.low[1,1:6], ci.u = tplot2.ci.high[1,1:6], xlim = c(0,0.5), horiz = TRUE, cex.names = 1.2, cex = 1.2, col = cols2[1,], las = 2, axisnames = FALSE)
mtext("Condition III", side = 3, line = 0, cex = 1.2)
mtext("Partisan cue", side = 3, line = 2, cex = 1.4)

barplot2(tplot1[2,1:6], plot.ci = TRUE, ci.l = tplot1.ci.low[2,1:6], ci.u = tplot1.ci.high[2,1:6], xlim = c(0,0.5), horiz = TRUE, cex.names = 1.2, cex = 1.2, col = cols1[1,], las = 2, axisnames = TRUE)
mtext("Condition II", side = 3, line = 0, cex = 1.2)
mtext("Deliberative cue", side = 2, line = 16, cex = 1.4)

barplot2(tplot2[2,1:6], plot.ci = TRUE, ci.l = tplot2.ci.low[2,1:6], ci.u = tplot2.ci.high[2,1:6], xlim = c(0,0.5), horiz = TRUE, cex.names = 1.2, cex = 1.2, col = cols2[1,], las = 2, axisnames = FALSE)
mtext("Condition IV", side = 3, line = 0, cex = 1.2)

#-------------------------- Table 2 --------------------------#

T2.upper <- tplot1.w.se[c(5, 4, 3, 2, 1, 6), ]

T2.lower <- tplot2.w.se[c(5, 4, 3, 2, 1, 6), ]

# print Table 2

xtable(T2.upper, digits = 1)
xtable(T2.lower, digits = 1)

#-------------------------- Figure 3 --------------------------#

point.ests <- cbind(tplot1[3, ], tplot2[3, ]) * 100

tplot1.ci <- cbind(tplot1.ci.low[3, ], tplot1.ci.high[3, ]) * 100

tplot2.ci <- cbind(tplot2.ci.low[3, ], tplot2.ci.high[3, ]) * 100

enames <-  rownames(point.ests)

par(mfrow =  c(1, 2)) 
par(oma =  c(3, 14, 0, 0), mar = c(0, 0, 3, 1)) 

plotCI(y = 1:6, x = rev(point.ests[,1]), err = "x", ui = rev(tplot1.ci[, 2]), li = rev(tplot1.ci[, 1]), axes = FALSE, ann = FALSE, ylim = c(0.7, 6.1), xlim = c(-30, 30), pch = 20, pt.bg = par(cex = 1.25), sfrac = 0.015, lwd = 1)
axis(side = 1, cex.axis = 1)
axis(side = 2, las = 1, at = c(1:6), cex.axis =  1, labels = rev(enames), las = 2)
mtext("No partisan cue", at = 0, line = 1, cex = 1.4)
abline(v = 0, col = "grey63")

plotCI(y = 1:6, x = rev(point.ests[,2]), err = "x", ui = rev(tplot2.ci[, 2]), li = rev(tplot2.ci[, 1]), axes = FALSE, ann = FALSE, ylim = c(0.7, 6.1), xlim = c(-30, 30), pch = 20, pt.bg = par(cex = 1.25), sfrac = 0.015, lwd = 1)
axis(side = 1, cex.axis = 1)
mtext("Partisan cue", at = 0, line = 1, cex = 1.4)
abline(v = 0, col = "grey63")

#-------------------------- Table 3 --------------------------#

dlexpdata$dem <- ifelse(dlexpdata$pid3 == 1, 1, 0)
dlexpdata$rep <- ifelse(dlexpdata$pid3 == 3, 1, 0)

# OLS

modl.ols <- lm(agree.num ~ video + rep + dem + video:rep + video:dem + party.cues:rep + party.cues:dem + party.cues + female + age + hispanic  +income + educ + video:party.cues + video:party.cues:dem + video:party.cues:rep, data = dlexpdata)

coefs.ols <- summary(modl.ols)$coefficients[,1:2]

# ordered logit

modl.ologit <- zelig(as.factor(agree.num) ~ video + rep + dem + video:rep + video:dem + party.cues:rep + party.cues:dem + party.cues + female + age + hispanic  +income + educ + video:party.cues + video:party.cues:dem + video:party.cues:rep, model = "ologit", data = dlexpdata)

coefs.ologit <- cbind(modl.ologit$get_coef()[[1]], sqrt(diag(modl.ologit$get_vcov()[[1]]))[1:length(modl.ologit$get_coef()[[1]])])

# make Table 3

Table.coefs <- matrix(c(t(rbind(coefs.ols[-1,], coefs.ologit))), ncol = 2, byrow = F)

rownames.inter <- c(t(matrix(rep(rownames(coefs.ols)[-1], 2), ncol = 2)))

rownames.inter <- paste(rownames.inter, c("coef","s.e."))

rownames(Table.coefs) <- rownames.inter

T3 <- Table.coefs[c(1:2, 19:22,3:8,27:28,23:26,31:32,29:30), ]

# print Table 3

xtable(T3, digits = 2)

#-------------------------- Figure 4 --------------------------#

set.seed(3000)

effects.by.pid <- list(NULL, NULL)

for (j in 1:2) { for (i in 1:3) {
  
  if (i == 1) {d <- 1; r <- 0} else{ if (i == 2) {d <- 0; r <- 0} else {d <- 0; r <- 1}}
  
  x.low <- setx(modl.ologit, video = 0, party.cues = (j - 1), rep = r, dem = d, female = 1, hispanic = 0, fn = list(numeric = mean))
  x.high <- setx(modl.ologit, video = 1, party.cues = (j - 1), rep = r, dem = d, female = 1, hispanic = 0, fn = list(numeric = mean))
  s.out <- sim(modl.ologit, x = x.low, x1 = x.high)
  
  change.approve.mean <- apply(s.out$`.->sim.out`$x1$fd[[1]], 2, mean)[5]
  change.approve.ql <- apply(s.out$`.->sim.out`$x1$fd[[1]], 2, quantile, p = 0.025)[5]
  change.approve.qu <- apply(s.out$`.->sim.out`$x1$fd[[1]], 2, quantile, p = 0.975)[5]
  
  effects.by.pid[[j]] <- rbind(effects.by.pid[[j]], c(change.approve.mean, change.approve.ql, change.approve.qu) * 100)
  
}}

effects.by.pid

title <- c("No partisan cue", "Partisan cue")

par(mfrow=c(1,2)) 
par(oma=c(1,1,0,0),mar=c(3,4,4,2)) 

newx <- 1:3

for (i in 1:2) {
  
  preds <- effects.by.pid[[i]]
  
  # plot
  plot(newx, preds[,1], type = 'n', ylim = c(-5, 20), bty = "n", xaxt='n', yaxt='n', ann=FALSE)
  axis(1, at = c(1, 2, 3), labels = c("Democrat", "Independent", "Republican"), cex.axis = 1.2)
  axis(2, at = seq(-20, 20, by = 5), cex.axis = 1.2)
  # add fill
  polygon(c(rev(newx), newx), c(rev(preds[ ,3]), preds[ ,2]), col = 'grey80', border = NA)
  abline(h = 0)
  lines(newx, preds[ ,1],  col = 'black')
  lines(newx, preds[ ,3], lty = 'dashed', col = 'grey12')
  lines(newx, preds[ ,2], lty = 'dashed', col = 'grey12')
  mtext(title[i], side = 3, line = 1, cex = 1.4)
  mtext("Treatment effect", side = 2, line = 3, cex = 1.2)
  abline(h = 20, lty = 3, col = "gray12")
  abline(h = 15, lty = 3, col = "gray12")
  abline(h = 10, lty = 3, col = "gray12")
  abline(h = 5, lty = 3, col = "gray12")
  abline(h = -5, lty = 3, col = "gray12")
}
